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It is thought that planets form from solid particles in a flattened, rotating, 99% gaseous nebula. 
These grains gradually coagulate into millimeter-to-meter sized aggregates which settle toward the 
midplane of the nebula. It is widely believed that the resulting dense layer eventually becomes 
gravitationally unstable and collapses into “planetesimals.” A new numerical model is presented to 
simulate the predominant processes (gravitation, vertical convection, and shear-driven turbulence) 
during the stage while the particulate material is still dispersed about the midplane of the nebula. In 
our previous work, particles were assumed to be spheres of a single radius; in the present work, 
particles are spheres of different radii. Results indicate that neither a broad nor a narrow distribution 
of particle sizes is likely to become gravitationally unstable. © 1995 American Institute of Physics . 


I. INTRODUCTION 

Dense regions of interstellar clouds collapse under their 
own gravity into young stars, surrounded by flattened disks 
of dusty gas. 1 Our own Solar system presumably originated 
from such a protoplanetary nebula. Subsequent collisions of 
asteroid-like “planetesimals” (followed by hydrodynamic 
capture of gas in the outer Solar nebula) led to their accretion 
into planets. 

The solid component of the nebula initially consisted of 
microscopic grains well-mixed with the gas. These grains 
would gradually coagulate into millimeter-to-meter-sized ag- 
gregates, which then settle into a thin layer at the central 
plane of the disk due to the vertical component of the Sun’s 
gravity. 2,3 It is widely believed 4 5 that as its density increases, 
this layer eventually grows gravitationally unstable and col- 
lapses directly into planetesimals. The criterion for gravita- 
tional instability may be written as 

Pp >M/(irR 3 ), (1) 

where p p is the bulk density of solid particles, M is the mass 
of the Sun, and R is the distance from the Sun. 6 

In recent years, however, it has become clear that plan- 
etesimals must form by other means, because settling of 
solid particles is a self-limiting process. 2 3,7,8 A radial pres- 
sure gradient causes the gaseous disk to rotate at slightly less 
than the Kepler orbital velocity. In contrast, the particle-rich 
layer is not supported by pressure, and so trends toward 
Keplerian motion. The resulting shear between the particle 
layer and the surrounding gas generates turbulence that stirs 
particles back into the nebula and reduces their bulk density 
below the stability limit. 

In order to simulate this situation, we have developed a 
computational fluid dynamics (CFD) model for multiphase 
flow in which the protoplanetary nebula is treated as a mix- 
ture of gas molecules and solid particles. The temporal and 
spatial evolution of the mixture is described by the 


Reynolds-averaged Navier-Stokes equations for each phase. 
The gas eddy viscosity is obtained from a Prandtl turbulence 
model adapted to the nebula flows of interest. The particle 
diffusivity is modeled by means of a Schmidt number ex- 
pressed as a function of particle size and density. 

Our initial results, presented in Cuzzi et ai? confirm 
that particles of a given mass and radius settle into a stable 
distribution about the midplane. In general, the lighter the 
particles, the thicker the layer. An interesting loophole re- 
mains, though. If particles of different sizes coexist, it is 
possible that the lighter particles may spread into such a 
thick layer that shear is reduced and turbulence is sup- 
pressed. Then the heavier particles could settle into a dense 
sublayer, which might become gravitationally unstable after 
all. 

The goal of the present work is to examine this possibil- 
ity by including particles of different sizes and compositions 
simultaneously. Our mathematical model is described in Sec. 
II, divided into six subsections detailing our assumptions, the 
model nebula, the governing equations, the perturbation 
technique, the Reynolds averaging method, and the turbu- 
lence model used to close the system. Section III describes 
the numerical technique, while Sec. IV presents the results. 

II. MATHEMATICAL MODEL 
A. Assumptions 

The assumptions used in modeling the protoplanetary 
nebula are outlined below. The gas phase is assumed to be a 
perfect gas, known to be predominantly hydrogen. The solid 
phases are each assumed to consist of spherical particles of a 
fixed mass and radius. Each phase is treated as a continuum. 
Momentum is exchanged between the gas and particle 
phases by means of the drag force. Each particle phase is 
independent of the other particle phases; they interact only 
through the intermediary of the gas phase. The settling of 
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particles is modeled by taking into account turbulence gen- 
erated in the gas by the presence of the particles. Collisions 
among particles are not taken into account, so the model is 
valid only for dilute suspensions; however, it contains no 
limitation on the mass loading ratio. The mean flow is as- 
sumed axisymmetric {3136= 0) and isothermal. Brownian ef- 
fects on the particle phase are neglected, as well as electro- 
magnetic forces. The self-gravity of the nebula is also 
neglected. Cuzzi et al . 8 demonstrate that this assumption is 
valid as long as condition (1) is not close to being satisfied. 


B. Nebula model 

The physical characteristics of the protoplanetary nebula 
are now introduced. A standard “minimum mass” circums- 
tellar nebula is assumed with total mass 0.0425 A/ 0 . The 
surface mass density oir) of the disk is of the form 

a(r) = cr 0 (r/AU) - 1-5 , (2) 

where AU represents 1 astronomical unit (1.49X10 13 cm). 
This leads to a surface mass density tr () ^ 1700 gem 2 at r = 1 
AU. 

The disk is also assumed vertically isothermal at tem- 
perature 

T(r) = 7' 0 (r/AU) -0 ' 5 , (3) 

where T 0 ^280 K. Because this is too warm for water ice to 
precipitate, the condensible mass fraction of 5.3 X 10 3 gives 
the solid particles a surface mass density of about 9 gem -2 
at 1 AU. 

Equations (2) and (3) above imply a gas density 

p g (r,z) = p tef (r/A\jy' U4 e\p[-(j/j G ) 2 ], (4) 

where z G is called the scale height. Note z G is always much 
greater than the depth of the particle layer, so the variation of 
gas density with height can be neglected. For example, 
£ G ^8X10 6 km at r=l AU in the baseline model. The gas 
density p ref at 1 AU is 1.4X10 -9 g cm -3 . The corresponding 
molecular mean free path is about 1 cm, so most of the 
particles are in the Stokes drag regime. 


C. Governing equations 

The equations describing the Solar nebula are expressed 
in a cylindrical coordinate system (r y 0 7 z) as 
Gas phase: 
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where t is time, p is density, and P is pressure. Here u, v, 
and w are respectively the components of velocity in the r 
(radial), 6 (circumferential), and z (vertical) directions. The 
subscripts g and p refer to the gas phase and particle phases, 
respectively. Here n is the number of distinct particle sizes; n 
is limited only by available computer resources; G is the 
gravitational constant, equal to 6.7X10 8 cm 3 g > S 2 . Thus 
the vertical and radial components of the Sun’s gravitational 
attraction become GMzIR 3 and GMrIR 3 , respectively. Fi- 
nally T f - j is the molecular stress tensor expressed in cylindri- 
cal coordinates, equal to r,y= fi(Uij + Uj /) where p,=0.001 
g/cm/s is the gas molecular viscosity. Equation (12) is pre- 
sented in conservative form for reasons given in Sec. III. 

The drag coefficient can be expressed as A p = (p g t p )~ l 
(Ref.8), where t p is the particle response time and depends 
on the ratio of the particle radius r p to the mean free path 
of a gas molecule. For (Epstein flow regime), 


l p 


Ps_ [p 
Pg c ' 


( 13 ) 


Here p s is the density of an individual solid particle, and c is 
the speed of sound in the gas. For r P>K (Stokes flow 
regime), 
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tp 3 P* C 0 |v p -v g | ’ 
where C D is the drag coefficient: 

C D — 24 Re ” 1 for Re^,<l, 

C d =24 Re ” 0 6 for l<Re p <800, 
C d = 0.44 for Re p >800. 

Here Re,, is the Reynolds number 
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where |v p — vj is the magnitude of the relative velocity be- 
tween the particle phase and the gas. 

D. Perturbation technique 

Because we are interested in relatively small variations 
of physical quantities over large distance ranges, it is appro- 
priate to solve the preceding set of equations in a perturbed 
form. The equations are expanded about a reference state as 
shown below. We seek a solution such that: u g = u [f 
v g = v 0 + v ]y w^w,, P*=Po+Pi, P = P 0 +P } , T=T 0 , 

Up=“p\i V p =V p0 + V pl , W p =w p y Pp=P p 1 . 

Substitution of this expansion in the exact equations 
leads to the following equations at zero order: 
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Here R g is the gas constant (8.31 43 X 10 7 erg/K/mole) and m 
is the mean molecular weight of the gas (2.34 amu). The 
quantity v K is the Keplerian orbital speed, and 77 is the frac- 
tional deviation of the gas from Keplerian rotation (propor- 
tional to the radial pressure gradient). To first order, the per- 
turbed state will satisfy the following equations in their exact 
form: 
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E. Reynolds averaging 

The perturbation quantities themselves are further ex- 
panded into means and variations, in the usual notation: 
p!=p + p\ u { = u + u\ V\ =v + v', W| == w + w'. By defini- 
tion, variation quantities vanish on the average: p' 

= 0, u' = 0, v* = 0, w' = 0, etc. This form is substi- 
tuted into the perturbation equations (23)-(30), which are 
then time averaged (Reynolds averaged) over an interval 
short compared to the evolution of the nebula but long com- 
pared to the time scales of the turbulence. 

Originally, we solved the vertical gas momentum equa- 
tion (8) for vv, the mean value of the gas vertical velocity 
perturbation. However, the Courant-Friedrichs-Lewy (CFL) 
condition required a very small time step for Eq. (8) to be 
numerically stable. We soon found that neglecting w elimi- 
nated the need for this CFL condition, and the time step 
could be increased an order of magnitude. Since w was al- 
ways small, setting vv=0 made no noticeable difference in 
the results. 

Since the gas density p g is nearly constant, its fluctuation 
will be very small, and any correlation involving p g is ne- 
glected. In the simulations presented in this paper, radial gra- 
dients (dldr) of perturbation quantities are also neglected by 
scaling. The resulting set of equations is shown below in the 
form in which they are solved. The overbar sign on all mean 
variables has been eliminated for clarity; 


d a 

T t p + ~r [(po+p)h-]=o, 
dt oz 


(31) 


du du p dP {) 2vi i(] + u‘ T eg 

dt dz Po(Po + p) dr r (Po + P) r 

1 a „ ^ (32) 


+ 


, , x, A p p p (u-u p ), 

(po + p) dz - P = l 


dv dv dv o du 0 u(Vq + v) 

- + W—--U — W — 

dt dz dr dz r 


+ 


* r6 


1 d 


■ + 


(p 0 +p)r (po+p) dz z 


- 2 A p p p (v 0 +v-v p() -v p ), 
p = i 


w = 0. 


(33) 

(34) 


Jt Pp + -^ { Pi’ w P + Pp w 'p ) = 0 ' 


(35) 


du 

~dt 


P . x ' m P 

+w »t 


du„ 2v p U p0 +Hp | ^ 


1 d 


■ + 


P p r Pn dz pzr 


p„w_ du n 

— T A p {p 0 +p)(u p -u), 

Pp dz 1 


(36) 


dv p dv p 

^T + Wp ^T~ 

dv p q 

U p dr Wp 


Up(v p() + Vp ) 


dv p0 


P p r p p dz p ~° 


_ Pp w p dVpp J_ d_ 

P p dz p p dz yi 'i”'n-r 

+ p)(Vpi) + Vp-V 0 -V), 


( v pP'p w 'p)~A p (p o 


(37) 


y t (Pp w p) + — (p p w p ) 

GM dr’ d 


. d 

Pp z+ ~dl 


~ 2 V Z (P P W P W P ) 


A p (po + p)p p {w -w) 


-A p (p i) +p)(p' p Wp-ppW' !i ). 


(38) 


In the viscous stresses, the superscript t stands for tur- 
bulent, while T stands for th e tota l of the molecular and 
turbulent stresses: r\ } = — P%u\u rjj = r t) + r^. The 
gas stresses are expressed as 
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where k is the kinetic energy of the gas turbulence. The 
particle stresses are 
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F. Turbulence modeling 

To close the system of equations, we adopted the formu- 
lation of Cuzzi et a/., 8 summarized as follows. Applying the 
Prandtl mixing length model to nebula flows, we express the 
turbulent eddy viscosity fi t as the square of a mixing length 
times the local gas density and vorticity: 

Here we take the mixing length equal to a constant c >,=0.045 
times the boundary layer thickness <5^0.02 tjv K /fl K . 
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FIG. 1. Geometry of nebula flow simulations. 


Correlations between particle density and particle veloc- 
ity are specified according to the gradient diffusion hypoth- 
esis. Since the particle density gradient is dominated by its 
vertical component, this implies 


t t r\ r r n t f * 

Pp“p = 0 ' Pp v P = 0 ’ Pp w p = ~s 


Uj_ 

■p g dz 
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Here S c is the Schmidt number, given by 


S r =(l+S,)Vl+3^;p g /2*, (43) 

where 5,=!!^ is called the particle Stokes number, and 
t p =(A p p g )~ 1 is the time constant of the drag. 8 The turbulent 
kinetic energy k is obtained from the eddy viscosity fx t 
through the relation 





k 


(44) 


where R 0 is the critical Rossby number (empirically set to 
80) and C p is a dimensionless coefficient (set to 0.09; Ref. 

8 ). 

Correlations between particle density and gas velocity 
are taken approximately equal to the corresponding particle- 
particle correlations. Consequently, the drag term involving 
the difference ( p p w p — PpWg) cancels out of Eq. (38) for 
the vertical momentum of the particles. 


III. NUMERICAL TECHNIQUE 
A. Algorithm 

The multiphase flow code is an extension of the two- 
phase flow code, with all particle variables expressed as ar- 
rays to take into account the different sizes of particles. The 
code employs dimensional variables, expressed in CGS 
units, but this can easily be changed. In this code, the time- 
dependent Navier-Stokes equations (31)-(38) are solved us- 
ing a time-marching approach. Time-dependent methods are 
commonly used in Computational Fluid Dynamics applica- 
tions to obtain a steady-state numerical solution of a fluid 
flow. By their use, a boundary-value problem is transformed 



FIG. 2. Model results at 1 AU for particles 60 cm in radius after 1 3 years of 
simulated time. The upper panel (labeled “G”) displays the gas velocity 
components u,v,w along with the particle density profile p p , while the 
lower panel (labeled “P") shows the velocity components u p ,w p for the 
particles. The diffusion velocity w diff and terminal velocity w F are defined in 
the text. 


into an initial-boundary value problem with unknown initial 
data (see Yee, Ref. 9, p. 127), for which efficient algorithms 
can be used. 

Implicit algorithms to solve partial differential equations 
with stiff source terms, like those in this paper, are still at the 
development stage (see Yee, Ref. 9, p. 118). Because our 
goal is to obtain a numerical solution to an as-yet unsolved 
problem, we use a well-developed explicit algorithm. We 
selected the 1969 MacCormack scheme, 10 because it has 
widely been used in aerodynamical simulations, although it 
has not been applied to astrophysical problems until this 
work. 7,8 We refer the reader to the work of Mendez-Nunez 
and Carroll, 11,12 who recently evaluated the MacCormack al- 
gorithm for atmospheric problems and demonstrated its ad- 
vantages with respect to numerical stability, diffusion, and 
dispersion. 

The radial and azimuthal momentum equations (32), 
(33), (36), and (37) are solved in nonconservative form. This 
form is employed instead of the more common conservative 
form in order to minimize roundoff errors in the particle 
velocities. In the conservative form, such errors may arise 
when the conservative variables p p u p ,p p v p ,p p w p are di- 
vided by the density p p to obtain the primitive variables 
Up,Vp,w p . With the nonconservative form, this is avoided 
since we solve directly for the primitive variables. 

However, we encountered a problem with the particle 
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FIG. 3. Similar to Fig. 2, but all particles are 10 cm in radius. Note change 
in scales. 

vertical momentum equation. Because of the terms contain- 
ing PpWp which originated from averaging advection terms, 
the conservative form (38) is actually more stable. In the 
one-dimensional case treated here, the particle continuity 
equation (35) and vertical momentum equation (38) are 
strongly coupled to each other, but coupled to the other two 
particle momentum equations (36) and (37) only through the 
value of the turbulent viscosity /jl { . Because of this situation, 
we have adopted a hybrid solution scheme. We solve (35) 
and (38) simultaneously in conservative form using first- 
order, upwind differencing for the advective terms (sacrific- 
ing some accuracy for the sake of stability), and second- 
order central differencing for the viscous terms. Then p p and 
w p are updated to the next time step, and the gas and particle 
radial and azimuthal momentum equations are solved in non- 
conservative form. 

Although this approach was stable, the algorithm tended 
to lose particle mass, due to excessive numerical diffusion 
arising from the first-order spatial differencing. Sophisticated 
numerical techniques that conserve particle (or species) mass 
to machine accuracy have been designed for unsteady mul- 
tidimensional aerosol problems, in particular by Toon et al . ] 3 
However, because we are looking for the steady-state solu- 
tion, a simpler approach was employed to conserve particle 
mass. 

We periodically applied an integrodifferential correction 
to compensate for the loss of particle mass. In a steady state, 
dpp/dt- 0, so that the particle continuity equation (35) and 
formula (42) give 



FIG. 4. Similar to Fig. 3, but for the 60 cm particles in a nebula containing 
equal masses of particles 10 and 60 cm in radius. Note change in scales. 


- - t t r - - dPp , 

Pp w p + P p ”p = c = Pp*p { 45) 

where C is a constant. Evaluating the left side of Eq. (45) 
above at z—0 gives C— 0, while integrating it leads to 

Pp( z ) = Pp(°) ex p( J o dzj. (46) 

In the numerical model, the value of p p ( 0) is determined 
from the total initial particle mass, and the profile of p p (z) is 
updated at each time step from Eq. (46) above. This tech- 
nique was found to be quite successful. 

B. Boundary conditions 

Our code was designed to treat flow in both the vertical 
and radial directions. Because so little is known about the 
particle settling, though, we have thus far restricted ourselves 
to one-dimensional calculations. Since the largest gradients 
lie in the vertical direction, we have neglected the complica- 
tion of radial flow for the time being. Accordingly, we define 
the numerical grid with either 102 or 202 rows in the vertical 
direction, but only three columns (IE =3) in the radial direc- 
tion. This geometry is depicted in Fig. 1. Zero-gradient con- 
ditions are imposed at both radial boundaries 7 = 1 and 1=3, 
so that all dependent variables are the same in each vertical 
column. Similar conditions are imposed on the horizontal 
rows K= 1 and K= 3, corresponding to symmetry across the 
midplane at K= 2. 
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FIG. 5. Similar to Fig. 4, but for the 10 cm particles in a nebula containing 
equal masses of particles 10 and 60 cm in radius. 



FIG. 6. Similar to Fig. 4, but after 51 years of simulated time. 


For the upper boundary conditions, we adapted the ap- 
proximate analytic solution of Nakagawa et al . 14 They were 
modeling a simple two-phase nebula, but we extended their 
results to multiple particle sizes. The results imply that as the 
particle density approaches zero far from the midplane, the 
density and velocity variations of the gas also vanish. In 
contrast, the particle velocity variations approach the limits 


u 


p 


S t 

l+sf ’ 


- vvk 

T+sJ’ 


w p =-zfl K S,. 


(47) 


We therefore imposed the above conditions on the density 
and velocity variations of each phase at the upper boundary 


K-KE. 


C. Code verification 

The nature of nebula flows does not permit direct com- 
parisons of simulations with experiment. Instead, computa- 
tions were verified by means of “internal” checks. First, a 
trivial test was performed. The new multiphase code was 
compared to our “old” two-phase flow code 8 for a mixture 
of gas with particles of a single size. Results for the two 
cases were identical. Next, these were compared to a multi- 
phase computation with all particles again of the same size, 
but now treated as two separate phases. The mass of each 
particle phase was half the mass of particles in the two-phase 
flow run. The results were still unchanged, as expected. 


IV. RESULTS 

Our first results are described below. The flow is a mix- 
ture of solid particles of individual density p s - 1.0 g/cm 3 
with nebular gas at r = 1 AU where the temperature is 280 K. 
Initially, the mesh consisted of 102 points with a vertical 
resolution of 600 km and a time step of 1/2 hour. After 11 
years of evolution, the mesh was refined by linear interpola- 
tion to 202 points with a vertical resolution of 300 km, and 
restarted with a time step of 1/4 hour. This resolution is still 
~ 100 times larger than the smallest eddies, and the time step 
is orders of 10 longer than their turnover time. In each case, 
the multiphase flow runs extend over about 13 years of simu- 
lated time, requiring about 2 hours of CPU time on a Cray 
Y-MP. Over the duration of the runs, residuals dropped by 
more than 4 orders of magnitude, indicating good conver- 
gence toward a steady-state solution. 

The results of these simulations are shown in Figs. 2-4. 
Figure 2 represents a two-phase run in which all particles 
were 60 cm in radius, while Fig. 3 displays a similar case 
where all particles were 10 cm in radius. These cases were 
selected because of the contrast in the thickness of the par- 
ticle layers. For comparison. Figs. 4 and 5 present a multi- 
phase flow containing equal masses of both species; Fig. 4 
refers to the particles of radius 60 cm while Fig. 5 applies to 
those of radius 10 cm. The total mass of all particles was the 
same in every run. 

The dot-dashed curves in each figure plot the bulk den- 
sity p p of the particle cloud; the scale is shown at the bottom 
of the upper panel. The solid and dashed curves in the lower 


Phys. Fluids, Vol. 7, No. 7, July 1995 


Champney, Dobrovolskis, and Cuzzi 


1709 








FIG. 7. Similar to Fig. 5, but after 51 years of simulated time. 


panel (labeled “P”) respectively show the radial and azi- 
muthal components u p and v p of the particle velocity, mea- 
sured with respect to the Keplerian speed v K (bottom scale). 
The corresponding gas velocities u and u are shown in the 
upper panel (labeled “G”), relative to pressure-supported 
circular motion i> 0 (top scale). 

The light dotted line in the lower panel represents the 
particle vertical velocity w p . The heavy dashed curve is the 
diffusion velocity 

— - — - — f-ij ^ d — 

defined as the effective velocity of the mass flow driven by 
the particle density gradient. 8 The good agreement between 
w diff and w p in the dense particle layer indicates that the 
particles are not settling appreciably. 

The heavy dotted curve plots the terminal velocity w F , 
the speed of fall for which the particle’s weight is balanced 
by drag: 

GM z - 

-jr Apo ~ tonZtpiz). (49) 

For the 60 cm particles of Figs. 2 and 4, w F differs from w p 
because particles falling from higher levels retain some mo- 
mentum. In contrast, w F is indistinguishable from w p for the 
10 cm particles of Figs. 3 and 5, indicating that the falling 
speed of the small particles continually adjusts to the local 
terminal velocity. 


Note that the particle density profile p p is much thicker 
in the 10 cm case (Fig. 3) than in the 60 cm run (Fig. 2). This 
is due to the greater diffusivity of the smaller particles. As a 
result, the profiles of azimuthal velocity v and v p are also 
thicker in Fig. 3 than in Fig. 2. Overall, the results of the 
multiphase calculations (Figs. 4 and 5) are fairly similar to 
the results of the two-phase runs (Figs. 2 and 3). As ex- 
pected, the velocity gradient in the multiphase case is inter- 
mediate between that in the two unimodal cases. This actu- 
ally causes the 10 cm particles in Fig. 5 to settle into a 
thinner layer than in Fig. 3. However, the distribution of 60 
cm particles in Fig. 4 is no denser than that in Fig. 2. This 
indicates that more realistic models including a broad distri- 
bution of particle sizes are not likely to be less stable in the 
Goldreich-Ward sense. 5 

In Fig. 4, we also observe numerical oscillations in the 
radial velocity u p of the 60 cm particles that were not present 
in the two-phase flow computations. We verified that these 
oscillations are reduced with a finer mesh, confirming that 
they are of numerical origin, and caused by our unwilling- 
ness to use an artificial viscosity. 

Note in each figure that drag between the particle-rich 
layer and the surrounding particle-poor nebula causes a ra- 
dial outflow of the gas. This outflow u peaks at —5,000 km 
above the midplane in Fig. 2, and at —30 000 km in Fig. 3, 
just above the respective particle-rich layers. In the multi- 
phase case, however, we observe two maxima in u corre- 
sponding to both particle sublayers. This effect could have 
significant implications for the radial transport of small par- 
ticles in the protoplanetary nebula and the compositional in- 
homogeneity of planetesimals. 

Finally, we extended the multiphase run from 1 3 years to 
51 years of simulated time. Lack of time and resources pre- 
vented our optimizing the time step, so we simply increased 
it by a factor of 10 and restarted the run. The results, dis- 
played in Figs. 6 and 7, are practically indistinguishable from 
those at 13 years plotted in Figs. 4 and 5. This demonstrates 
that a steady state has indeed been reached. 
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